function psi_t_solver(m,ω,xmin,xmax, Npoints,x0, p0,sigma, time)
#Basis
b_position = PositionBasis(xmin, xmax, Npoints)
b_momentum = MomentumBasis(b_position)
# operators
x = position(b_position)
p = momentum(b_momentum)
px = transform(b_momentum, b_position) # position to momentum
xp = dagger(px)
K = LazyProduct(xp, p^2/2m, px) #kinetic energy
V = ω^2*x^2*m/2 # potential energy
H = LazySum(K, V) #Hamiltonian
psi0 = gaussianstate(b_position, x0, p0, sigma) #assume initial state is a gaussian function
t, psi_t = timeevolution.schroedinger(time, psi0, H);
return t,psi_t
end
function SplitStep(ψ, V, p, t_points, m; potential_evolution=x_coeff, kinetic_evolution=p_coeff)
#=
ψ: initial wavefunction
V: time independent potential
p: momentum
T: total time
Δt: time-step
=#
Δt = t_points[2]-t_points[1]
output = Array{Complex}(undef, length(t_points), length(ψ))
# Creates the empty output array. 1st dimension is time, including initial wvfcn, 2nd is position.
for (i, t) in enumerate(t_points)
output[i,:]=ψ[:] # Inputs the current content of ψ into the output
ψ=fft(potential_evolution(V, Δt).*ψ)
ψ=ifft(kinetic_evolution(p, Δt, m).*ψ)
ψ=potential_evolution(V,Δt).*ψ
end
output
end
x_coeff(V::Vector{ComplexF64}, Δt::Float64) = exp.(-im*V*Δt/2); # Creates the coefficient for position evolution
p_coeff(p::Vector{ComplexF64}, Δt::Float64, mass::Float64=1.) = exp.(-im*(p.^2)*Δt/(2mass)); # Creates the coefficient for momentum evolution
function split_step_method(m,ω,x_min,x_max, Npoints ,x0, p0,sigma, time)
#=
ψ: initial wavefunction
V: time independent potential
p: momentum
T: total time
Δt: time-step
=#
x = Vector{ComplexF64}(LinRange(x_min,x_max,Npoints));
Δx = x[2]-x[1]
p = 2π.*Vector{ComplexF64}(fftfreq(length(x), length(x)/abs(x[1]-x[end])));
#p or po
ψ = (sqrt(Δx)/(sqrt(sigma)*(π)^(1/4)))*exp.(im*p0*(x.-0.5x0)-((x.-x0).^2)./(2sigma^2))
V = ω^2*x.^2*m/2;
output = SplitStep(ψ,V,p,time, m);
return output
end
